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Abstract 

Horizontal wells present several advantages compared to vertical wells. However, analytical 
models developed so far can only describe the pressure behavior in single-layer reservoirs, 
while there is no available solution for multilayer stratified reservoirs with horizontal wells. 
Based on Newman’s product and source functions, this work aims to achieve an analytical 
formulation for stratified reservoirs under single-phase flow. The developed model’s accuracy 
was verified by comparison to a numerical flow simulator. The proposed formulation was 
also used to estimate the reservoir equivalent permeability. 
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1. Introduction 

Interest in developing analytical solutions 
comes from the possibility of inferring 
several reservoir features from the pressure 
data collected during a transient test. 
Well testing is a procedure that aims to 
estimate reservoir properties by analyzing 
pressure transient data when the reservoir 
is submitted to single-phase flow. Initially, 
the reservoir underlying fluid is produced, 
leading to a pressure drop. This stage is 
denoted as drawdown. During the next 
stage, which is called buildup, the wellbore 
is shut. Pressure, then, increases. Current 
interpretation methods allow, for instance, 
the computation of the reservoir equivalent 
permeability and the identification of skin 


*1 am corresponding author 
Email addresses: 

renanvbl@aluno.puc-rio.br (Renan Vieira 
Bela), sinesio@puc-rio.br (Sinesio Pesco), 
abelardo.puc@gmail.com (Abelardo Barreto Jr.) 

Preprint submitted to Elsevier 


effect and the reservoir outer boundary 
[1, 2, 3], 

Well testing may be performed in both 
vertical and horizontal wellbores. The 
latter present some advantages, such as 
its superior performance in reservoirs with 
small thickness or high vertical permeability 
[4], Horizontal wells may also produce at 
higher flow-rates, which reduces operational 
costs. Besides, fewer horizontal wells are 
required to produce the same amount of oil, 
when compared to vertical wells [5]. 

Despite those gains, existing analytical 
models are only capable of describing the 
pressure behavior in single-layer reservoirs 
with horizontal wells [6, 7]. To the 
best of the authors’ knowledge, the only 
available solution for multilayer reservoirs 
with horizontal wells assumes that there is 
formation crossflow; that is, fluid transfer 
between layers may occur [8]. However, 
there is no analytical formulation for 

November 24, 2019 



systems here there is no communication 
between layers other than through the 
wellbore, which are denoted as commingled 
or stratified reservoirs. 

Thus, the main purpose of this work is 
to achieve an analytical model for pressure 
transient behavior in stratified systems 
under single-phase flow. The proposed 
model was based on the existing solutions 
for single-layer reservoirs. Furthermore, 
the reservoir equivalent permeability was 
estimated using the suggested formulation. 
The accuracy of the developed model was 
assessed by comparison with a commercial 
flow simulator. 

In the following sections, a brief recap 
regarding the previous achievements is 
conducted. Then, the considered reservoir 
model and the proposed formulation are 
presented. Next, the results are displayed, 
followed by the main conclusions of this 
work. 

2. Previous Achievements 

The first analytical model for single-layer 
reservoirs with horizontal wells was 
presented by Goode and Thambynayagam 
[6]. Their solution was achieved by 
successively applying Laplace and Fourier 
transforms into Darcy’s law. It accounts 
for anisotropy and assumes uniform influx, 
that is, flow-rate per unit length is constant 
along the well. However, this assumption 
does not match the physical reality of the 
wellbore, which is more accurately depicted 
by the infinite conductivity hypothesis 
(pressure inside the well is constant). 
To overcome this issue, they stated that 
pressure at an infinite conductivity wellbore 
may be approximated by the average 
pressure computed from a uniform influx 
based model. Moreover, buildup pressure 


is determined applying the superposition 
principle. 

The real space solution for single-layer 
systems was reached by Daviau et al. 
using source functions and Newman’s 
product [7]. The analytical model 
also assumes uniform influx. However, 
infinite conductivity is estimated from the 
equivalent pressure point technique. They 
identified two distinct radial flow regimes 
and briefly commented the influence of 
wellbore hydraulics and positioning. 

Following their work, anisotropy effects 
were included into the real space solution 
by Ozkan et al. [4]. Furthermore, 
a means for computing the reservoir 
permeability and skin is proposed based 
on pressure derivative analysis. Further 
details regarding early and late time radial 
flows are provided and the equivalent 
pressure point technique is endorsed as an 
alternative to estimate pressure behavior 
in an infinite conductivity wellbore using a 
uniform influx model. 

To obtain a more accurate pressure 
response for the infinite conductivity 
hypothesis, Rosa & Carvalho divided 
the wellbore into several uniform influx 
segments [9]. Pressure equality inside 
those segments induces a linear system that 
allows the determination of flow-rate in 
each segment. The linear system must be 
updated at every time step, since flow-rate 
in each segment varies in time. 

Up to four distinct flow regimes may 
be observed in reservoirs with horizontal 
wells under single-phase flow [3]. At early 
time, pressure behavior is similar to a radial 
flow in an infinite acting reservoir with a 
vertical well, where the reservoir thickness 
equals the horizontal wellbore length. After 
the reservoir vertical boundaries begin to 
influence the pressure response, a linear flow 
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develops, followed by a late radial flow. 
This second radial flow occurs when flow in 
the horizontal plane becomes more relevant 
than in the vertical axis. The last observed 
regime is the boundary governed flow, which 
starts after the measured pressure begins to 
reflect the reservoir lateral boundaries. 

The solution in Laplace domain is 
applicable for different vertical boundary 
conditions, such as no-flow boundary, gas 
cap or aquifer [10]. Pressure maintenance 
systems may be easily identified due to their 
characteristic pressure derivative signature. 
The model assumes once more uniform 
influx along the well. Thereby, pressure 
profile for infinite conductivity wells, may 
be obtained using either the equivalent 
point or the average pressure technique. 

A more realistic wellbore model is 
achieved by accounting for frictional and 
hydraulic losses inside it [11]. Then, 
pressure change is no longer a linear 
function of the production flow-rate. 
Similarly to the infinite conductivity case, 
flow-rate profile along the well changes in 
time, but stabilizes after a short transient 
period. 

The only analytical solution for 
multilayer reservoirs with horizontal 
wells assumes formation crossflow in all 
layers [8]. The wellbore is considered to 
penetrate only one reservoir layer and fluid 
transfer at layer interfaces is depicted by a 
combination of reflection and transmission 
coefficients. Thus, pressure response is 
progressively influenced by layers further 
from the wellbore, as the transient pressure 
pulse propagates throughout the reservoir. 

This overview evidences that there is 
currently no formulation for multilayer 
stratified reservoirs with horizontal wells. 
Thus, this work attempts to adapt existing 
analytical models for single-layer systems 



Figure 1: Considered Reservoir Model 


so that they might be extended to the 
multilayer case too. 

3. Proposed Formulation 

Figure 1 displays the considered reservoir 
model. It is assumed a multilayer stratified 
reservoir with a multilateral horizontal 
well, such that each layer is perforated 
by exactly one well ramification. Since 
no formation crossflow is considered, this 
wellbore configuration plays an essential 
role so that pressure response reflects 
properties from all layers. Otherwise, only 
layers penetrated by the wellbore would 
influence the measured pressure profile. 

Anisotropy effects are accounted for, 
as the vertical permeability is different 
than the horizontal permeability. In fact, 
the suggested formulation is more general 
and can also represent the case when 
permeabilities in x and y directions are 
not the same. Hereafter, all computations 
assume a consistent set of units. Other 
relevant simplifying hypothesis are listed 
below: 

• Homogeneous, anisotropic, laterally 
infinite reservoir 

• Reservoir initially in equilibrium 

• Only the horizontal wellbore sections 
are perforated 
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• Pressure is the same in all layers, apart 
from the hydrostatic effect 

• Slightly compressible fluids with 
constant viscosity 

• Constant production flow-rate 

• Gravitational and wellbore storage 
effects are neglected 
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Since the reservoir is vertically bounded, 
the images method must be applied to 
obtain the pressure response in z [7]: 


3.1. Solution for Single-Layer Reservoirs 

Applying Newman’s product and the 
sources/sinks method in a single-layer 
reservoir, the instant pressure drop owing to 
all source point productions along the well is 
obtained by the product between the three 
directional pressure changes and the source 
intensity [4, 7]: 


A P ina (x,y,z,t) = 

y^—AP x {x, t)AP y (y, t)AP z (z, t) ^ ^ 

Equation (1) is achieved assuming 
uniform influx. The production flow-rate 
is represented by q, while L stands for the 
wellbore length, cj) denotes the porosity 
and c t is the total compressibility. Setting 
the well heel at a given point (xo,yo,zo), 
and defining the hydraulic diffusivity as 

rj = -, directional pressure change in x 

4>yc t 

is given by [9]: 
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In equation (4), dz stands for the smallest 
distance between the wellbore and a vertical 
boundary. Pressure change at a given time 
is defined as the difference between the 
initial pressure and the measured pressure 
at this time. It is determined by integrating 
all instant pressure variations defined by 
equations (1) to (4) and including the skin 
effect: 
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In y direction, the contribution of all 
points along the wellbore must be accounted 
for [4]: 


In equation (5), S represents the 
mechanical skin factor, as defined by 
Hawkins [12]. 

3.2. Formulation for Multilayer Reservoirs 
The formulation for multilayer reservoirs 
starts by applying equation (5) into a given 
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layer j. Initially, it will be assumed that 
there is no formation damage, simply for a 
matter of convenience. Thus: 


APj(x,y,z,t) = —^— 

Lj^jCTj 

} ( 6 ) 

/ A P xj (x, T)AP yj (y, r)AP zj (z, r)dr 

o 

Equation (6) is achieved under the 
assumption that layer flow-rates remain 
constant during the test. If layer properties 
are remarkably distinct, flow-rates change 
in time. Then, the qj term cannot 

be written out of the time integral in 
equation (6). However, based on the 
formulation for multilayer reservoirs with 
vertical wellbores [1, 2], this rate transient 
period is not expected to last long. 
Directional pressure changes in layer j are 
determined analogously to equations (2) 
to (4). Thereby, layer flow-rate may be 
determined as: 


Qj 
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where the Oj term in equation (7) is defined 
as: 
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Thus, summing the flow-rates in all layers 
to obtain the total production flow-rate: 
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By model hypothesis, pressure is the 
same in all layers, except for the hydrostatic 
effect. Hence, equation (9) may be 
rearranged yielding: 


A P(x,y,z,t) = q 
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Last, the damaged region in each layer 
must be accounted for. In order to do that, 
the reservoir overall skin factor is given by 
[1]: 
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Besides, the equivalent permeability in 
the zx vertical plane must be defined: 
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In equation (12), L T stands 
wellbore total length. Thus: 
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Equation (13) represents the pressure at 
a given point for a reservoir subject to 
the assumptions mentioned above. No 
restriction is made regarding the well 
positioning. That is, the developed 
formulation accounts for well off-centering. 
Properties such as permeability, porosity, 
skin factor and well length can also be 
different for each layer. The strongest 
hypotheses made to achieve equation (13) 
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consist of assuming uniform influx and 
constant layer flow-rates. 

Nevertheless, the proposed formulation 
can easily represent an infinite conductivity 
wellbore using either the average pressure 
technique [6] or the equivalent pressure 
point [4, 7]. In this work, the former method 
was used, with some adjustments due to the 
considered multilateral wellbore model. 

Each well ramification was segmented at 
100 equidistant points. This number of 
wellbore segments was arbitrarily chosen. 
Then, for a given time, 100 pressure 
values were computed setting x = z = 
r w . In each layer, the a 3 term was 
evaluated applying, in y direction, the 
wellbore length corresponding to those 
100 segments. Wellbore pressure was 
determined by taking the average between 
those 100 computed pressures. Figure 2 
represents an example where each wellbore 
ramification was divided into 4 different 
segments. The generalized expression for 
the average pressure in a wellbore divided 
into N s distinct segments is given by: 
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Figure 2: Schematics showing the points used to 
compute the wellbore average pressure 


The constant flow-rates hypothesis, in its 
turn, is not expected to result in significant 
errors, since it is quickly satisfied [1]. 
During buildup, the superposition principle 
is used to determine the pressure response 
[6, 8]. Buildup pressure change is defined 
as the difference between the pressures 
measured at the last drawdown point and 
a given buildup time. 

3.3. Pressure Derivative Behavior During 
Each Flow Regimes 

In single-layer reservoirs with horizontal 
wells submitted to single-phase flow, up to 
four distinct time regimes may be observed: 
early radial, intermediate linear, late radial 
and boundary governed [3]. 

The early radial flow occurs while the 
pressure transient pulse has not reached 
the reservoir vertical boundaries and flow 
beyond the well tips is negligible. It 
resembles the flow through a vertical well in 
an infinite acting reservoir whose thickness 
equals the horizontal well length [10]. 
During this flow regime, pressure derivative 
with respect to the natural logarithm 
of time (hereafter referred to simply as 
pressure derivative) reflects the geometric 
mean of the permeabilities in x and z 
directions [3]. Assuming once again a 
consistent set of units, pressure derivative 
in multilayer stratified reservoirs may be 
approximated as: 

tWt K w (15) 

5 In (t) L T {k zx ) eq 

Since the considered reservoir model 
accounts for distinct layer properties, 
the pressure transient pulse propagates 
differently in each layer. Therefore, the 
intermediate linear flow might not develop 
simultaneously in all layers. Hence, the 
transition period between early and late 
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Case 

ho 

(cP) 

Q 

(m 3 /d) 

Layer 

0 

h 

(m) 

L 

(m) 

dz 

(m) 

k H 

(mD) 

k z 

(mD) 

S 

A 

3,6 

1100 

1 

0,18 

15 

150 

5 

500 

170 

2,3 

2 

0,25 

20 

150 

10 

750 

190 

3,3 




1 


10 

100 

5 

310 

100 

0,9 

B 

1,0 

800 

2 

0,27 

10 

120 

5 

520 

170 

2,3 




3 


10 

150 

5 

740 

250 

3,3 


Table 1: Reservoir Properties 


radial flow is more accurately depicted by 
a 3D transitional flow rather than a linear 
flow regime per say. For this reason, 
pressure derivative presents no a priori 
characteristic signature. 

After some time, flow in the vertical 
plane becomes negligible and the late 
radial flow regime is reached if the 
reservoir is sufficiently wide. During this 
period, pressure transient front propagates 
radially throughout the horizontal plane 
[10]. Pressure derivative, then, reflects the 
permeabilities in directions x and y: 

^- (16) 

<91n(t) h T (k xy ) eq 

where: 


y ] hj yjk x j k 


yj 


(K)e, = 3 - i — 7 - ( 17 ) 

iLt 

The last flow regime starts when the 
reservoir outer lateral boundaries are felt 
by the wellbore. In this work, a laterally 
infinite reservoir is considered. Thus, the 
boundary governed flow regime is never 
achieved, and the last identifiable flow 
regime is the late radial. 


4. Validating the Analytical Model 

To verify the accuracy of the analytical 
solution detailed in the previous section, 


results obtained for two cases were 
compared with a finite difference-based 
commercial flow simulator. In both cases, 
the reservoir was modeled using a cartesian 
grid with 5.0 km wide (x direction) and 2.5 
km long (y direction). These dimensions 
are sufficiently large for the laterally infinite 
reservoir hypothesis to hold. Oil model 
used was blackoil. The employed numerical 
grid was progressively refined so that the 
highest discretization takes place closest 
to the wellbore, while regions closer to 
the reservoir lateral boundaries are more 
roughly represented. 

In both cases, it was considered that 
well drilling and completion change the 
permeability around it to 100 mD in all 
directions, representing the skin effect. 
Grid blocks containing the wellbore were 
built with 0.6 m width and 0.6 m 
height. Permeability inside those blocks 
was changed to 100 mD to represent the 
damaged region. In the analytical solution, 
the skin zone was modeled defining the skin 




Figure 3: Damaged Zone Representation 
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Figure 4: Pressure and Pressure Derivative Profile for Case A 


radius required so that the cross-sectional 
damaged region would present the same 
area as the numerical simulator, yielding 
a skin radius of 0.35 m. Figure 3 
illustrates this procedure. In both cases, 
wellbore radius was set as 0.11m and total 
compressibility as 1.72xl0~ 4 (kgf/cm 2 )" 1 . 
Other relevant reservoir properties may be 
found in Table 1, which considers the 
Brazilian field units set. 

Besides the pressure response, pressure 
derivative with respect to the natural 
logarithm of time was also analyzed for both 
numerical and analytical data. Derivative 
profile was computed using Bourdet 3-point 
algorithm [13]. During buildup, derivative 
was taken with respect to Agarwal’s 
equivalent time [14], 

4-1. Results and Discussion 

Pressure and pressure derivative profiles 
for cases A and B are displayed in figures 4 
and 5. The theoretical pressure derivative 
levels associated with the early and late 
radial flow regimes (equations (15) and 
(16)) are also plotted in dashed lines. On 
the other hand, the boundary governed 


regime is not observed during the test. 
This is an immediate consequence of the 
laterally infinite reservoir assumption in the 
analytical model and the large reservoir 
dimensions employed at the numerical 
simulation. 

For both cases, a close agreement 
between analytical and numerical data is 
identified. This suggests that the constant 
flow-rates hypothesis required to develop 
the proposed formulation was indeed 
quickly reached and the early rate transient 
period did not reflect in significant errors, 
as expected. Moreover, the good observed 
match also indicates that the adjustments 
proposed to compute the pressure response 
in multilateral infinite conductivity wells 
using the average pressure method were 
successful. 

The main divergences in both cases 
were noticed at early times. However, 
pressure derivative profile suggests that the 
analytical formulation is, in fact, more 
accurate at early times than the simulator, 
as the former attains the theoretical early 
time pressure derivative level more precisely 
than the latter. 
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Figure 5: Pressure and Pressure Derivative Profile for Case B 


This fact may be explained by some 
issues regarding the numerical simulation. 
Pressure is computed by the simulator 
at the center of each block. Therefore, 
simulated data are sensitive to the employed 
grid refinement [15]. Besides, the damaged 
region in the simulation grid is modeled 
by cubes, although the skin zone is in 
fact cylindrical. Furthermore, the simulator 
represents the wellbore as a line source. 
These reasons might explain the observed 
mismatch at early times. 

The same effects are more intensely 
observed during buildup, resulting in more 
relevant divergences between analytical 
and numerical data at early buildup 
times. Moreover, since layer permeabilities 
are different, there might occur an 
inside-wellbore crossflow between layers 
at early buildup times. This effect is not 
accounted for by the analytical formulation, 
and is another possible explanation for the 
early buildup divergences. 

In case B, wellbore lengths are different in 
each layer. Thus, early buildup divergences 
might also be related to the fact that 
the zero-rate pulse does not simultaneously 


reach the well toe in all ramifications. 
This inside-wellbore rate transient possibly 
justifies the more significant errors during 
early buildup times in this case. 

Moreover, the early radial flow is shorter 
than in case A. This is explained by 
the smaller layer thicknesses and wellbore 
lengths in case B. Early radial flow lasts 
while flow beyond the well tips is negligible 
[3] or the reservoir vertical boundaries 
are felt by the well [10]. Thus, the 
transitional flow starts later in case A, 
where layer thicknesses and wellbore lengths 
are larger. The same reservoir and wellbore 
features also explain why the late radial flow 
develops sooner in case B. 

4-2. Determining the Reservoir Equivalent 
Permeability 

Based on the two constant pressure 
derivative levels, two distinct estimates for 
the reservoir equivalent permeability may 
be determined. Those estimates, which are 
associated with the early and late linear 
flow regimes, are computed using equations 
(15) and (16), respectively. The computed 
results for cases A and B may be seen in 
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Case 


( kzx 

)eq 


{J^xy 

)eq 

Real 

Est. 

Error (%) 

Real 

Est. 

Error (%) 

A 

335 

350 

4,7 

643 

647 

0,7 

B 

318 

330 

3,6 

523 

525 

0,3 


Table 2: Estimated Equivalent Permeabilities 


Tabic 2. All permeabilities were computed 
with excellent accuracy. 

The higher errors presented by the 
computed permeabilities in zx plane are 
related to the fact that pressure derivative 
attained values slightly lower than the 
theoretical level during early radial flow. 
Moreover, it is important to mention that, 
in a real field test, early radial flow may not 
be identified in reservoirs where the storage 
effect is relevant. Regardless, permeability 
in the horizontal plane could be determined 
with great precision in both cases. 

5. Conclusions 

An analytical formulation for stratified 
reservoirs with horizontal wells under 
single-phase flow was developed. The 
proposed model derives from the existing 
solution for single-layer systems and is 
based on Newman’s product and source 
functions. To achieve the suggested 
formulation, it was assumed uniform 
influx along the well and constant layer 
flow-rates during the test. Moreover, the 
average pressure technique was adapted to 
estimate the pressure profile at an infinite 
conductivity multilateral wellbore based on 
a uniform influx model. 

The suggested formulation was validated 
after comparison to a commercial numerical 
flow simulator. Results indicates that the 
rate transient period was too short to 
imply in significant errors. Furthermore, 
the proposed adjustments at the average 


pressure method proved to be efficient for 
the considered wellbore model. 

Estimates for the reservoir equivalent 
permeability were obtained from the 
proposed formulation. For each case, 
two distinct permeabilities may be 
computed, one representing the equivalent 
permeability in the vertical plane 
perpendicular to the well axis and another 
related to the properties at the horizontal 
plane. Both estimates presented low errors, 
showing that the developed solution may 
be helpful in providing information about 
the reservoir. 


Nomenclature 


Ct = Total compressibility 

dz = Smallest distance between the well 

and a vertical boundary 

h = Thickness 

k = Permeability 

L = Wellbore length 

P = Pressure 

q = Flow-rate 

r w = Wellbore radius 

S = Mechanical skin factor 

t = Time 

rj = ffydraulic diffusivity 
// = Viscosity 
cj) = Porosity 

r = Silent time integration variable 
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Subscripts 

eq = Equivalent 
ins = Instant 
j = Related to layer j 
x = Related to x direction 
y = Related to y direction 
z = Related to z direction 
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